Oscillations in two-species models: tying the stochastic and deterministic approaches 
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We analyze general two-species stochastic models, of the kind generally used for the study of pop- 
ulation dynamics. We show that the conditions for the stochastic (microscopic) model to display 
approximate sustained oscillatory behavior are governed by the parameters of the corresponding de- 
terministic (macroscopic) model. We provide a quantitative criterion for the quality of the stochastic 
oscillation, using a dimensionless parameter that depends only on the deterministic model. When 
this parameter is small, the oscillations are clear, and the frequencies of the stochastic and deter- 
ministic oscillations are close, for all stochastic models compatible with the same deterministic one. 
On the other hand, when it is large, the oscillations cannot be distinguished from a noise. 
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It is well known that the dynamics of many systems 
of two species can display an oscillatory behavior in the 
populations of both agents. This happens in predator- 
prey systems in models of measles epidemics Q, in 
chemical systems such as those exemplified by the Brusse- 
lator Q , etc. These systems are usually modelled by a set 
of two coupled ordinary differential equations, which are 
assumed to represent a macroscopic level of description 
of the system. Oscillations can appear in these mod- 
els as limit-cycle solutions to the equations. However, 
it frequently happens that the macroscopic model only 
has damped oscillatory solutions, even though the mod- 
elled system displays sustained oscillations in its popula- 
tions in the same region of parameter values. Examples 
of this are not uncommon in population dynamics (see, 
e.g., the discussion in with regard to predator- prey 
and measles problems) . It has often been noted that the 
stochastic counterpart of these models — assumed to rep- 
resent a more microscopic level description of the same 
system — usually do display a kind of sustained oscilla- 
tory behavior, with a frequency very similar to the one of 
the damped solutions of the differential equations @,[1,0| 
(see Fig. Q] for an example based on a susceptible- infected 
epidemic model). These oscillations are said to be gen- 
erated by the demographic, or intrinsic, noise The 
problem is that stochasticity precludes a clear-cut defi- 
nition of "oscillations" for such systems. Therefore, the 
comparison between the results of the stochastic and de- 
terministic approaches is often made on a qualitative ba- 
sis. 

In this Letter we address the problem of sustained os- 
cillations in stochastic models, in an attempt to charac- 
terize the oscillatory regime. We show that the conditions 
for well defined oscillations are given by the parameters 
of the corresponding macroscopic (deterministic) model 
only, disregarding the details of the microscopic (stochas- 
tic) one. In other words, this general result proves that 
the conditions are the same for any stochastic model that 
corresponds to the same macroscopic one. To this end, 
we define a criterion of "quality" of oscillatory behavior 



for a large class of stochastic two-species systems and we 
show, by means of a van Kampen expansion of the master 
equation, that the information given by the determinis- 
tic system — embodied in a deterministic parameter — is 
enough to provide good bounds on this quality. In other 
words, we show that the quality of the oscillations is 
only weakly dependent on the details of the demographic 
noise. Moreover, it is shown that oscillations become 
clear if and only if the deterministic parameter vanishes. 
We also suggest a heuristic value for the quality below 
which one can be almost certain that the evolution of 
both populations does "look" oscillatory. 

We consider systems of two populations, A and B, de- 
scribed by stochastic variables m{t) and n(t). The state 
of the system is defined by the joint probability P(m, n; t) 




FIG. 1: Deterministic dynamics (smooth lines) and one 
stochastic realization (fluctuating lines) of an SI epidemic 
model (susceptible and infected, respectively A and B). The 
dynamics includes birth and death processes in both popu- 
lations, and contagion. A self-limiting intraspecific compe- 
tition mechanism is implemented as in [|, E3, with a to- 
tal system size O = 5 x 10 5 . Transition rates are: Tio = 
2{bA(j>A + bBA 4>b){1 — <I>a — 4>b), T_io = dA<j>A, To-l = dB4>B, 
T_ii = 2 P <f> A <t>B (see Eq. (TlJ). 
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that the system has m individuals of species A, and n in- 
dividuals of species B. The transition from a state with 
(m, n) individuals to a state with (m+i, n+j) individuals 
takes place at a rate: 



T(m + i,n + j\m,n) = f(n)T i:j ( — , -), 



(1) 



where — k < i < k and —k < j < k. is a scale pa- 
rameter that governs the fluctuations of the stochastic 
evolution. Its precise definition depends on the system, 
but one chooses it in such a way that for large fl the 
fluctuations are small. It usually represents the volume 
containing the reactants in chemical systems or the 
available resources in biological ones [8|. The constant 
k gives the maximal number of elements that can ap- 
pear, or disappear, from a given population at each step 
of the dynamics. The most common choice are one-step 
processes, with k = 1. 

The evolution of the probability P(m, n; t) is given by 
the master equation 
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where, as in the rest of this Letter, the summation indices 
run from — k to k. 

Except for a few simple cases, this equation is ex- 
tremely difficult to solve exactly. For this reason many 
methods have been devised to look for approximate solu- 
tions. Perhaps the best known, and most applied, is the 
van Kampen expansion 0]. In the following we sketch 
the main steps leading to the series solution (a detailed 
account can be found in van Kampen's book [9(). 

If one assumes that, at time zero, the system is in a 
state where both populations have well defined macro- 
scopical values, P{m, n) — 6(m — mo)8(n — rto), with the 
initial values of order O(fi), it is reasonable to expect that 
at later times P(m, n) will have a sharp peak at some po- 
sition of order 0(fl) (in both populations), and a width 
of order 0(p}/ 2 ). That is, the fluctuating populations 
will satisfy m = fl(f)A + \/^<£a and n = fl(j)B + V&£,b, 
where the variables (j> represent the "macroscopic" evo- 
lution, while the stochastic variables £ represent fluctu- 
ations around them. Replacing this in Eq. ([2]), equating 
terms of the same order in il and adequately rescaling 
the time, one obtains, for the leading order: 



<t>A = 22iTjj((t>A,<t>B) = C A {(j>A,(j>B)i 
ij 

4>B = y^JTjj(<f>A,<i>B) = Cb^a^b)- (3) 

ij 

These equations, called deterministic or macroscopic, are 
usually the starting point of many models of chemical and 
biological systems. They are generally written down from 



macroscopic considerations of the population dynamics, 
disregarding its individual level origin. To analyze the 
differences between the stochastic (individual level) and 
the deterministic (population level) approaches one usu- 
ally chooses a stochastic model that gives the right deter- 
ministic equations. In the limit of infinite size (f2 — * oo), 
Eqs. ([2]) are also satisfied by the average populations. 

The deterministic equilibria are obtained by solving 
the system Ca(4>a, 4>b) — Cs(0yi,0B) = 0, and their 
stability is studied by means of a linear stability analy- 
sis. When the system is close to a stable equilibrium, its 
evolution can be approximated by that of a damped os- 
cillator. In the underdamped regime, the damping factor 
and the frequency of oscillation are, respectively: 
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with 



Ae/2, 
A(l- 

e=\T\/VA, 



74), 



(4) 



(5) 



where A is the Jacobian of C, and T its trace: 

A = Ca,aCb,b — Ca,bCb,a, 

T = C a ,a + Cb,b, (6) 

where C^j = The underdamped regime is there- 

fore given by the condition e < 2. We show below that 
this parameter, which depends only on the parameters of 
the macroscopic Eq. ([3]), plays a fundamental role in the 
characterization of the oscillations of stochastic origin. 
Notice also that the number of oscillations observed in 
the characteristic time 7 _1 depends only on e (for small 
e, it is just 2/e). 

The following order in the van Kampen expansion gives 
the evolution of 11(^,^5,*), the joint probability func- 
tion of the fluctuations, in the form of a Fokkcr-Planck 
equation. To look for oscillations in the fluctuations it 
is easier to work with the equivalent Langevin equations, 
as shown by McKane [8|: 

£,a — ~Ca,a£,a — Ca,b£,b + LA{t) 

£b = — Cb,aCa — Cb,b£b + L B {t) (7) 

where L-A(t) and Ls(t) are delta-correlated Gaussian 
noises of zero mean, satisfying (LA(t)LA(t')) = DA${t — 
t'), (L B (t)L B (t')) = D B S{t-t r ), and {L A {t)L B {t')) = 
DABS(t — t'). The noise intensities are given by: 



Da{4>a,4>b) = 2_^i Tij{4>A,4>B), 
Duty A , 4b) = ^j 2 T y (0A>s), (8) 
Dab(4>a, 4>b) 



/JjTij[<t>A,<t>B)- 



By Fourier transforming Eqs. ([7]) it is straightforward 
to obtain the power spectrum of the fluctuations around 
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FIG. 2: u>\ and Q% for the same SI model as in Fig. [T] The 
points correspond to a uniform scanning of a portion of phase 
space: b A = b B A = 0.1, d A G (0,0.2), d B G (0,0.5), p G 
(2d B ,3.2). The full lines show the bounds of Eq. JT2J), while 
the dashed one corresponds to u>d- The arrows point to the 
values corresponding to the parameters used in Fig. [3] 



the deterministic equilibrium [8| . In the following we con- 
centrate on population A. The corresponding expressions 
for population B are obtained by exchanging A and B in 
all the subindices. The average power spectrum of (a is 



F A +LU 2 



{l-Co 2 ) 2 +Co 2 e 2 ' 



(9) 



where 



Cj 2 = uj 2 /A, 

@A B^B + C% B^A — %Ca,bCb,bDab 



F a = 



AD, 



(10) 



We stress that Fa (and correspondingly Fg), through its 
dependence on Ca , Cb , Da and Db , depends ultimately 
on the transition probabilities that define the model. 

It is straightforward to see that (Sa(w)) is either mono- 
tonically decreasing or it has a single maximum at 



u? A = -F A + V(Fa + I) 2 - e*F A . 



(11) 



The condition of positivity for the argument of the square 
root gives the region in phase space where the power spec- 
trum has a single maximum. Notice that for e < \[2 this 
condition is fulfilled regardless of the exact dependence of 
Fa on the parameters of the model. It can also be proved 
that <jj 2 A satisfies [121 ]: 



1 - e 2 /2 < Cj\ < 1 



(12) 



(these bounds seem to be tight). In particular, this im- 
plies that in all the possible stochastic models that lead 
to the same deterministic equations (same C's, different 
D's) the position of the maximum can only vary within 
a finite range, that shrinks with e. 



For small e, Cj 2 a tends to 1, which means that not only 
the frequencies of possible oscillations for both popula- 
tions become close, but also that they become close to 
u>d, the frequency of the damped oscillations of the de- 
terministic model (which also tends to 1 as e — > 0). It 
is in this regime that the populations show the coherent 
dynamics characteristic of stochastic oscillations. This 
motion will be further characterized below by the qual- 
ity of the spectrum peak. Figure [2] shows Cj\ and u>% as 
functions of e for the SI model presented in Fig. [IJ for 
a wide range of system parameters. The bounds given 
by Eq. (| 1 2|) are shown by continuous lines. Each point 
represents the normalized squared frequency for one set 
of parameters, for both populations. The deterministic 
frequency u)^ is also shown, to emphasize the difference 
between the three frequencies present in the system. 

When V2 < e < 2 there can be some stochastic models 
for which no peak is present in S(uja) or S(u>b)- And, 
for some values of Fa or Fb, it can happen that the 
power spectrum of any population has a maximum even 
if e > 2, i.e. even when the deterministic system does not 
display damped oscillations (see Fig. O all the points to 
the right of e — 2 correspond to systems with a peak in 
the spectrum of the susceptible (A) population, no peak 
in the infected (B) one, and no damped oscillations in 
the deterministic model). These two features show that 
the peaks of the stochastic power spectrum on the one 
hand, and the deterministic damped oscillations on the 
other, are not necessarily closely related. 

The above discussion establishes the conditions for the 
existence of a peak in the power spectrum of one or both 
populations. That is, for the existence of a preferred 
frequency in their dynamics. But, should all peaks in 
the power spectrum be regarded as "oscillations"? The 
answer to this question is certainly negative, and leads 
one to look for a criterion to quantify how close a time 
series is to an oscillatory movement. This can be done by 
defining the "quality" of the oscillation as a measure of 
the sharpness of the peak. We propose one such measure 
in the following. 

Given a power spectrum of the form ([9]) we define the 
quality of a peak at Lu pe ak as 

U peak (S A(.Up ea k)) 



Qa(w) 



J (S A (oj))doj 



(13) 



This quantity is dimensionless and scale invariant. It 
is related to Fisher's kappa, which measures the non- 
stationarity of a signal, given its periodogram ll|. For 
functions with only one peak, Qa increases as the peaks 
grows. For power spectra of the form (|9|), Qa can be read- 
ily calculated (using that / (SA(w))duj = (£1), see ['J, ;: 



}a{oja) = 



LU A e 



(Ld A -l) 2 +u 2 A e 2 \F A + 1 



F A 



(14) 



The quality Q A diverges as e vanishes, regardless of the 
exact dependence of F A on the parameters of the model. 
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FIG. 3: Two stochastic realizations of the SI model, with dif- 
ferent qualities. The insets show the corresponding analytical 
average power spectra. Only the infected population is shown. 
The arrows in Fig. [2] point to the corresponding frequencies: 
uia ~ ujb ~ i^>d for the good quality case shown in (a), and 
uia ^> lob for the bad quality one shown in (b). Q. = 10 J . 



Therefore, one can assure that the corresponding time 
series will look oscillatory when e is sufficiently small (see 
Fig. [3] for an example of this) . In such a case, we have 
already shown that the frequencies of both populations 
are very close, and also very close to the frequency of the 
deterministic damped oscillations. 

Could it also happen that, for large values of e, when 
the frequencies of populations A and B can be rather 
different, one gets very sharp peaks? It can be shown 
that this is not the case by giving bounds of Qa that 
depend solely on e [121 ] : 



1 
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2 /4 



< Q A (u) < — • 



(15) 



The upper bound shows that, when e is not small, the 
peak cannot be arbitrarily sharp. On the other hand, 
the lower bound shows that, when e is small the peak is 
sharp for all the stochastic counterparts of a deterministic 
model. In Fig. 2] we illustrate this by showing several 
values of Qa and Qb for the SI model, along with the 
corresponding bounds. 

One practical question remains: what is the critical 
quality value above which one can be sure that the time 
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FIG. 4: Qa and Qa as functions of e for the SI model of 
Fig. [T] The points correspond to the same portion of phase 
space as in Fig. [2] The lines show the upper and lower bounds 
of Eq. (IS). 



series will indeed "look" oscillatory? As it is to be ex- 
pected, the continuous nature of Q precludes a conclusive 
answer. From exhaustive observations of different mod- 
els we find that when Q(uj pea k) > 1, the oscillations are 
well defined and notably different from a noisy evolution 
(see Fig. g]) . 

In summary, we have shown that, by defining a quality 
measure, one can quantify the "oscillatory look" of a time 
series. Interestingly, we find that oscillations are present 
only when e is small. This means that, given a deter- 
ministic model, one can know, using the bounds (|15[) . 
whether the time series given by any stochastic counter- 
part of the model will look oscillatory or not. In addition, 
we have shown that, when oscillations are clear, the cor- 
responding frequencies of both populations will be close 
to each other and to the frequency of the damped oscil- 
lation of the deterministic system. 

Given that our conclusions are based on the analysis 
of the first two terms of the systematic van Kampen's 
expansion of the master equation, they are exact only in 
the limit f2 — > oo. These analytical results, nevertheless, 
compare well with the numerical observations made on 
finite systems. More details about the validity of the 
expansion for finite systems will be given elsewhere 12j. 
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